arXiv: 1505.00797v2 [gr-qc] 13 Nov 2015 


Nonlinear interactions between black holes and 
Proca fields 


Miguel Zilhacd, Helvi Witek^, Vitor Cardoso*’* 

t Departament de Ffsica Fonamental & Institut de Ciencies del Cosmos, Universitat 
de Barcelona, Marti' i Franques 1, E-08028 Barcelona, Spain 
** Department of Applied Mathematics and Theoretical Physics, Centre for 
Mathematical Sciences, University of Cambridge, Wilberforce Road, Cambridge CB3 
OWA, UK 

" School of Mathematical Sciences, University of Nottingham, University Park, 
Nottingham NG7 2RD, UK 

* CENTRA, Departamento de Ffsica, Instituto Superior Tecnico, Universidade de 
Lisboa, Avenida Rovisco Pais 1, 1049 Lisboa, Portugal 

* Perimeter Institute for Theoretical Physics Waterloo, Ontario N2J 2W9, Canada 

E-mail: mzilhao@ffn.ub.es, h.witek@damtp.cam.ac.uk, 
vitor.cardoso@ist.utl.pt 

Abstract. Physics beyond the Standard Model is an important candidate for dark 
matter, and an interesting testing ground for strong-field gravity: the equivalence 
principle “forces” all forms of matter to fall in the same way, and it is therefore natural 
to look for imprints of these fields in regions with strong gravitational fields, such as 
compact stars or black holes. Here we study General Relativity minimally coupled to a 
massive vector field, and how black holes in this theory lose “hair”. Our results indicate 
that black holes can sustain Proca field condensates for extremely long time-scales. 


1. Introduction 

Black holes (BHs) are strongly gravitating objects that solve the equations of General 
Relativity (GR) and which are predicted to arise as the end-state of gravitational 
collapse. There is now overwhelming evidence that the universe is populated by millions 
of BHs, in mass ranges of a few solar masses to gigantic, billion solar-mass objects. 
Supermassive BHs are tightly connected to the evolution of their host galaxy and may 
prompt or quench star formation in the entire galaxy through complex interactions 
between the BH and its environment. 

Black holes play a crucial role in fundamental physics, as the “elementary particle” 
of gravity. This is largely due to the fact that both in GR and in many other 
modified theories of gravity, BHs are surprisingly simple objects characterized by only 
few parameters. In a series of rigidity, uniqueness and no-hair theorems M, it has 
been established that stationary BHs in Einstein-Maxwell theory (in four spacetime 
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dimensions) belong to the Kerr-Newman family of solutions and, thus, are completely 
characterized by their mass, angular momentum and (electric) charge. 

Uniqueness properties together with the demonstration of mode stability kicked off 
the success story of BHs as viable objects in the universe. While rigorous mathematical 
proofs of BH stability against generic perturbations are still underway [I], it has been 
shown that subextremal BHs are linearly stable against a wide class of perturbations 
involving massless fields IMS These perturbative results are now being complemented 
by studies concerning the nonlinear stability of BHs [11] [12], a property that is crucial 
for the understanding of highly dynamical scenarios. 

Black holes can be employed as sensitive “scales” to investigate ultra-light 
fundamental fields predicted in beyond-standard model physics [IMS]. The latter has 
become possible with the realization that the paradigm “BHs in four-dimensional GR 
are stable” is no longer true in the presence of massive bosonic fields. Instead, BHs may 
suffer from a superradiant or “BH-bomb” type instability mm, as has recently been 
proven rigourously [18] [§] This instability requires two ingredients: (i) Superradiance, 
the amplification of low-frequency waves when they scatter off spinning BHs P3|- For 
monochromatic waves of frequency ujr scattering off a spinning BH with angular velocity 
at the horizon, the condition is 

0 < lur < mQ H , ( 1 . 1 ) 


where m is the azimuthal “quantum” number of the wave. This provides a classical 
mechanism to amplify the impinging field at the expense of the BH’s rotational energy, 
(ii) Confinement of the field in the vicinity of the BH. As a consequence of these 
two ingredients, the field is successively amplified leading to an exponential growth 
of amplitude. 

Massive bosonic fields with mass parameter /i = ms/h create such a confining 
cavity, thus yielding superradiant instabilities |22}j24]. The timescale of the process is 
regulated by the coupling 


%Mn = 7.45 • 10 9 
ch 


M 

Mo 


rriR 


T// 2 i , ( 1 - 2 ) 

eV/c 1 J 

between the BH mass M and the field’s mass i%. This coupling is too large for 
astrophysical BHs (with masses ranging between a few solar masses up to 1O 1O M 0 ) 
and standard model particles - consider, for example, the Higgs boson for which 
M/r ~ 10 21 ... 10 3 ° - and yields instability timescales that are much longer than the 
age of the universe {24]; however, the “BH-bomb” mechanism becomes relevant if 


f Note, that extremal BH solutions are unstable [9], but it is unlikely that these solutions are ever 
attained in realistic astrophysical environments m- 

§ At the onset of the instability a branch of new, hairy BH solutions has been found for complex 
scalars. These configurations evade standard no-hair theorems and essentially interpolate between Kerr 
BHs and boson stars HSU20]- The time-development of the superradiant instability in an adiabatic 
approximation, and the appearance of scalar condensates was investigated both for real and complex 
scalars in m■ 
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Mji ~ 0(1) (in natural units). Thus, astrophysical BHs are sensitive to ultra-light 
bosonic helds in the mass range 10~ 23 eV < mgc 2 < lCU 10 eU. 

These fields appear naturally when more fundamental theories that attempt to 
marry gravity with quantum physics are compactified to four spacetime dimensions. 
Depending on their interpretation, these fields may play a role as modifications of GR, 
the simplest of which are scalar-tensor theories, where the scalar field may acquire an 
effective mass due to interactions with the environment [25:, 2D.- Similarly, there are 
extensions of GR involving additional vector fields such as Einstein-Aether theory ra. 
Horava-Lifshitz gravity [28], or higher dimensional gravity when compatified to four 
spacetime dimensions [2H|, to name but a few explicit examples. Although these 
additional vector fields are typically massless, environmental effects may induce an 
effective mass term as described above for scalar-tensor theories. Understanding the 
properties and phase-space of BHs in alternative models of gravity is an active field of 
research. For example, no-hair theorems in generic scalar-tensor theories have been 
discussed in Sotirious’ contribution to this Focus Issue ra- Of direct interest for 
the applications in the present manuscript is Bekenstein’s series of papers [31, [32] in 
which it was shown that static and stationary BHs cannot be bestowed with massive 
scalar, vector or spin-2 fields. This implies that, e.g., the Schwarzschild solution is 
the unique static solution of the Einstein-Proca system. An exception are the recently 
constructed “hairy BHs” that appear at the onset of the instability and evade these no¬ 
hair theorems [19] [20]. In the context of particle physics, the “axiverse” scenario suggests 
the existence of an entire landscape of axion-like particles, i.e., ultra-light pseudo-scalar 
fields that become massive through spontaneous symmetry breaking or nonperturbative 
effects [T3HT5] . Another natural extension of the standard model involves additional 
vector fields. In particular, compactifications of string theories (or their low-energy 
incarnations) yield a hidden U( 1) gauge group whose symmetry may be broken via 
a Higgs-like mechanism PS [33]. This spontaneous symmetry breaking may endow 
the hidden “photon” with a mass, in a manner similar to the standard-model Higgs 
mechanism that is responsible for the mass of the W ± and Z bosons. 

In summary, superradiant instabilities can turn astrophysical BHs into 
“laboratories” to investigate the properties of these fundamental fields. Potentially 
observable signatures include (i) gaps in the Regge plane, i.e. the BHs’ mass-spin 
phase-space [13], 2T]; (ii) modifications of the inspiral dynamics and characteristic, 
monochromatic gravitational wave (GW) emission [34l [37] : or (iii) gravitational 
radiation induced by a bosenova collapse [38]. 

With this motivation in mind, the phenomenology of massive fields in BH 
backgrounds has been an active field of research HZl Unfortunately, the timescales 
involved are so large that most studies have focused on perturbative calculations of 
the instability. Recently, we started a long-term programme aimed at understanding 
the nonlinear development of fundamental fields around spinning BHs, the first step of 
which focused on massive scalars [31]. A natural next step—and our purpose herein—is 
the extension of these studies to massive vector, i.e. Proca fields. Frequency-domain, 
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linearized calculations in the background of Schwarzschild BHs [391 ®! and slowly- 
rotating BH spacetimes mi m revealed a hydrogen-like frequency spectrum in the 
small-coupling limit. These results have been used to put new stringent bounds on the 
allowed mass of the hidden “photon”, p 7 < lCT 20 eId (42J. The growth rates for the 
superradiant instability of massive vector fields are up to four orders of magnitudes 
larger than those of the massive scalar BH and linearized, time-domain studies have 
predicted timescales of r = 3 • 10 3 ^|4 ~ O.OIstf- in the best case scenario 031 • For this 
reason Proca fields represent an appealing model to explore the “BH-bomb” instability 
at the fully nonlinear level. 

The present paper is the first installment of a series attending this challenging 
issue. Here, we focus on developing the formalism and presenting a full-blown, 3 + 1- 
dimensional numerical code capable of evolving BH spacetimes coupled to massive 
vector fields. Our results have been obtained through two different, independent 
implementations. For the sake of testing and benchmarking our codes we restrict 
ourselves to the evolution of massive vector fields in a nonrotating spacetime. The 
evolution of spinning geometries will be considered elsewhere. 

The paper is organized as follows: In section [2] we introduce our model and present 
its time evolution formulation. We derive constraint-preserving initial configurations 
in section [3] In section [4] we provide our formalism to extract information about the 
emitted radiation. In section [5] we give a brief summary of our implementation and 
present the numerical results. We finish with our conclusions in section [6} Throughout 
the paper we use natural units, i.e., c = 1 = G and H — 1. Greek letters /i, u,... refer 
to spacetime indices running from 0,..., 3, whereas Latin letters i, j,... denote spatial 
indices running from 1,..., 3. 

2. The model — coupling Proca fields to gravity 

2.1. Action and equations of motion 

We consider an action coupling a vector field X p with mass my = /iyh to gravity dans] 



( 2 . 1 ) 


where W pu = X p X u — V u X p . Note, that we consider solely the gravitational coupling to 
massive vector fields, because we ultimately aim at exploring superradiant scattering, 
a purely gravitational effect that is independent of the kinematic mixing between the 
hidden sector and visible photons. The resulting equations of motion (EoMs) are 


X v W pu + liyX 11 = 0 , 

Rfiu 2 9fj.i>R = 21 'V/j.p 


( 2 . 2 ) 

\g^w pa w pa + 4 (2X P X U - g^X p X p )( 2.3) 


g pv R = 2 W pp W u p 


In contrast to Einstein-Maxwell theory (where fly — 0), the vector X p is not a pure 
gauge field but plays the role of a fundamental, physically meaningful field. This implies 
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that the Lorenz condition 
V'.f,=0 ( 


(2.4) 


has to be satisfied, as can be seen from ( 2 . 2 ) 


2.2. Formulation as a Cauchy problem 

To explore the nonlinear interactions between Proca fields and BHs we have to track 


their dynamics numerically. Therefore, we recast the EoMs (2.2) and (2.3) as a time 


evolution problem. We begin by foliating the spacetime into 3-dimensional spatial slices 
E t with metric 7 ^. The 3-metric defines the projection operator 


Yv = 


(2.5) 


where nh denotes the t im elike unit vector (normal to E t ) with normalization = — 1 . 
The full, 4-dimensional spacetime metric can then be expressed as 

ds 2 = g^dx^dx 1 ' = — (a 2 — /3 l /3i) d t 2 + 2 r y i j/3 l dtdx j + 7 ^ da ; 1 dad , ( 2 . 6 ) 

where the lapse function a and shift vector /T describe the coordinate degrees of freedom. 
To complete the characterization of the full spacetime we define the extrinsic curvature 

Kij = — ^ Yt ~ £p) lij • (2-7) 

We now split the Proca field into its scalar X^ and 3-vector A) potential given by 

Xi = Yi X /x , and X^ = -n^X ^. ( 2 . 8 ) 


x j> x ji + ^(7/ ‘Tfp 


with 


In analogy to Maxwell’s theory we introduce an “electric” and “magnetic” field, 

Ei = YiW /iV n v , Bi = Yi *W lu ,n v = Y k D :] X k (2.9) 

where E^n 11 = 0 and = 0 hold by definition. As we will see later, the former 

definitions provide a time evolution equation for the 3-vector potential, while we employ 
the latter purely as short-hand notation. The W /IU tensor is reconstructed from 


fiu ^fiEy n v E^ 


D^X V — D u X fJ/ , 


( 2 . 10 ) 


where D M is the covariant derivative with respect to the 3-metric. 

The definition of the extrinsic curvature and the W /IU tensor prescribe the time 
evolution of the 3-metric 7 ^ and (spatial) vector potential X t , whereas the dynamics of 
the model are determined by the EoMs ( |2.2 ) and (2.3). Their spatial projection yields 
the evolution equation for the “electric” field and extrinsic curvature. Additionally, the 


Lorenz condition (2.4) prescribes the evolution of the scalar potential X^. 

Proca’s equations (just like its Maxwellian counterparts) give rise to constraints 
which the evolution equations are subjected to—the familiar “Gauss” constraint. 
However, analysing the evolution partial differential equations (PDEs) reveals the 
existence of a mode with zero characteristic speed. In other words, a constraint violation 
(which is always present in numerical simulations) will not be propagated away, but will 
rather stay in the numerical domain. It is therefore desirable to modify the Proca 
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evolution equations such that all characteristic speeds are non-zero and the constraints 
are darnpe <1 Gnided by numerical simulations of BH binaries in Einstein-Maxwell 
theory we have introduced an additional variable Z and damping parameter k > 0, 
which will help stabilising the numerical evolution [471 [49] . We recover the original 


system (2.2), (2.3) and (2.4) by setting Z — 0. In summary, the GR-Proca model 


evolves according to 

dtlij = -2aK i:j + (2.11) 

dtXi = — a (77,; + D % X^) — X^D^y. + /GW;, (2-12) 

d t E { = a (KE i + D*Z + v ^X i + e ijk DjB k ) - e ijk BjD k a + CpE\ (2.13) 
d t Kij = — DiDjOi + a (Rij — 2K ik K k j + /7/v i? j + 


2 a 


EiEj - ]pijE k E k + B t Bj - ]p ij B k B k - /lyXiXj ) ,(2.14) 


dtX* = - X'Dta + a (KX^ - D t X l - Z) + CpX h (2.15) 

dtZ = ol ( DiE 1 + fiyXfj) — kZ ^ + C,pZ , (2.16) 

where C denotes the Lie derivative. Additionally, the evolution is subject to a set of 
constraints given by 

U = R - K.jK'i + K 2 - 2 (E l E t + B^, + ^ (X 2 + A*A))) = 0 ,(2.17) 
Mi = MK tJ - D t K - 2 (e ljk E'>B k + ^X^) = 0 . (2.18) 


In practice, we adopt a free evolution scheme, in which the constraints (2.17), (2.18) 


are only solved for the initial data which will then be evolved. When setting up the 
relevant initial data, we also need to satisfy the “Gauss” constraint 

£ = D i E i + tiX (j) = 0, (2.19) 


which is enforced during the evolution through (2.16) 


Throughout the numerical simulation we will not actually use the ADM-York type 


formulation, given by (2.11)—(2.16), which is an ill-posed Cauchy problem, but we will 
employ the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) scheme [5(1, 5Tj together 
with the moving puncture gauge [52l f54j . This reformulation involves a conformal 


decomposition together with a constraint addition rendering it into a well-posed initial 
value problem [55!, 56]. The explicit expressions are rather lengthy and therefore given 
Appendix A 


m 


2.3. Flat space limit and Yukawa suppression 

Minkowski geometry is a good approximation at large distance, and provides a useful 
benchmark for our results, in particular for the time and spatial decay of the fields. The 
ansatz 




Qidr, / Q2df,0,0 


( 2 . 20 ) 


A similar procedure for the GR evolution sector also exists and is known as the “Z4” formulation (H 
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where Qi = Qi(t,r), produces the equation of motion for T = Qi — Q 2 , 

T" - T + - ^ r2 + = 0 , (2.21) 

ry* 

where ' and ' denote derivatives with respect to t and r. We find the harmonic solutions 

* = e ~ lut (CiJ- 2 fir^J/4-uj 2 ^ + C 2 Y_ 2 (iryjfi-utyj , (2.22) 

where Y, J are Bessel functions of the first kind and C % are constants. 

Another flat space solution (in 3 + 1 form) is given by 

E r (r ) = Ce~ kr (4 + - 1 cos ut (2.23) 

\r 2 r J 

X r (r) = -C^e~ kr (\ + -'\ sin ut (2.24) 

/iy \r 2 r J 

XA,(r)—C—^ -coscat, (2.25) 

hv r 

where cu = ^///y — fc 2 . 

Both solutions teach us that at large distances the fields have a Yukawa-like 
dependence on the radial distance for low-energy modes (cu < fi y) and an oscillatory 
behaviour for high-energy modes. Guided by the flat-space solutions we expect to 
recover a Yukawa-like dependence for quasi-bound states around BHs, because by 
definition these modes are trapped in the vicinity of the BH. 

2-4- Summary of perturbative results 

The above results concerning Yukawa suppression in Minkowski backgrounds are in fact 
shown to arise also in a perturbative framework of (low-amplitude) Proca fields around 
BHs. These results will be important as benchmarks for our numerical implementation 
of the nonlinear BH-Proca system. For this purpose, we briefly summarize the most 
important results on the subject. 

In contrast to Einstein-Maxwell theory, massive spin-1 fields propagate three 
physical degree of freedom. These come in the shape of an axial mode (labelled with 
S — 0) and two polar modes with spin S = ±1. Note, that the Proca field also 
propagates a monopole, i.e., spherically symmetric mode corresponding to S — +1. 

In a perturbative approach, the field equations are linearized in a fixed-BH 
background and any function (e.g., the vector potential) can then be Fourier-decomposed 
~ e~ lU}t . When boundary conditions are imposed, the held equations become an 
eigenvalue problem for the characteristic frequencies, which become complex quantities, 
u = ojr + noi- Frequency domain analysis of the Proca equation in the background 
of Schwarzschild BHs [39, jUH 15TH59] and slowly rotating BH spacetimes [H] [42j have 
revealed two classes of solutions that depend on the particular choice of the boundary 
conditions at spatial infinity. One type of solution corresponds to quasi-normal modes 
(QNMs) [6], while the second class are quasi-bound states (QBSs), which are particular 
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of massive fluctuations and, essentially, decay exponentially at large distances. Quasi¬ 
bound states require confined states, as we discussed previously, and only occur for 
Ur < 

For small mass couplings ilf/iy -C 1, the QBS are well described by hydrogen-type 
spectra mi. 


Ur ~ /iy 


(Mfl y) 2 \ 

2 (l + n + S + l) 2 ) ’ 


(2.26) 


where n > 0 is an overtone number and l is the angular quantum number [391141 j . The 
growth or decay rate of the QBSs is encoded in the imaginary part of their frequencies 
approximated by the analytic expression 


oji ~ ( ma/M — 2r + /iy) (M py) 41+5+2S , 


(2.27) 


in the slow-rotation regime [3Tj. In the limit a/M —>■ 0 it reduces to the results found 
in reference [39] for the Schwarzschild case. The superradiant instability is strongest in 
nearly extremal backgrounds and numerical computations in the time domain found a 
maximum growth rate of Mu>i ~ 5-1CT 4 (for dimensionless spin parameter a/M = 0.99) 
which corresponds to an e-folding timescale of r ~ O.Ols/^ [ 33] . 

The backscattering off the spacetime curvature gives rise to oscillatory power-law 
tails whose frequency is determined by the mass of the held 


0 ~ t p sin(/i V t), (2.28) 

with the exponent p = — (l + 3/2 + S) at intermediate late times M/r v ( 

and a universal decay p = —5/6 at very late times •C tp y [33] [60], 61]. Thus, 

for monopole and dipole fluctuations the held is expected to decay, respectively, as 
sin(/ry t) and t _ ( 5 / 2+fir ) sin(/xy t) at intermediate times and to follow the universal 
behaviour t~ 5//6 sin(/iyt) at very late-times, if the “tail-stage” is not superposed by the 
exponential, but slowly decaying bound-state stage. Generic initial data will excite all 
possible modes. Whether tails or QBS dominate the signal depends on the type of initial 
data and at which stage the signal is being observed. 


3. Initial data 

In this section we describe the construction of initial configurations representing a BH 
in the presence of a massive vector held. Although we only evolve nonrotating BH 
spacetimes in this paper, our construction in principle allows for rotating solutions. We 
intend to explore helds in the nonperturbative regime. Then, an important ingredient 
for long-term stable nonlinear simulations is the construction of constraint-satisfying 
initial data ]34j . For the purpose of solving the initial data problem we employ the 
York-Lichnerowicz conformal decomposition [62, S3] 

= + Wy-Vf, E i = -4,-<-E i , x i , = i,- e x 4 ,. 


(3.1) 
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The constraints become 


n = Aib--R+ -'—AvAa - —K 2 + 

8 8 -07 y 19 


^ r ,2 , /4 


12 


4-0 7 


ij + x t x 3 


+ (L* - D,X,) , 

A4j = ZUA^ - + 2& (tbjXi - DiX^j - 2 l j 2 r X <fj X l , 

,2 


£ — DiE 1 + fiyX^ , 


(3.2) 

(3.3) 

(3.4) 


where the hatted operators A and D % are constructed from the conformal metric 7 ^. 
These PDEs are in general difficult to solve, so we now simplify the constraint equations 
using the following assumptions 


% = Sij , K = 0 , Xi = 0, E { = -S ij djV. (3.5) 

yielding 

0 = D 3 A\ , (3.6) 

0 = M + w AiiAii + W s>ia>vdiV + ’ {3J) 

0 = AV - i^Xt . (3.8) 


In the following we will outline two different procedures to solve these equations. 


3.1. Simple Initial Data solutions 


A straight-forward way to solve the constraints (3.6)—(3.8) is provided by further 
imposing X$ = 0. Equation (3.6) has the known Bowen-York solution [64], whereas 
the Laplace equation AV = 0 is solved by a potential V ~ With these solutions 


the constraint equations are now decoupled, and one can solve (3.7) numerically for the 
conformal factor ijj using standard methods. 

If we additionally assume time-symmetric data, i.e. A t] = 0, the momentum 
constraint (3.6) is satisfied trivially and we find exact, analytic solutions for 
equations (3.8) and (3.7) (remember we are imposing X^ = 0) given by 


rz (1 M \ 2 C 2 _ C 

~ \ + 2 r) 4r 2 ’ V ~ r 1 


(3.9) 


where C is a constant. We term this a “Simple Initial Data” (SID) solution. Using our 
ansatz (3.5) we can read off the initial profile for the “electric” field, 


E i = il>~ 6 C^ , (3.10) 

where x l = (x,y,z). This construction is fully analytic, thus very simple to implement 
and run. In the Einstein-Maxwell case, /iy = 0, this reduces to an exact (static) 
solution describing a single Reissner-Nordstrom BH with charge C. In the present 
case, however, C has to be interpreted as the amplitude of the Proca field instead of 
an electromagnetic charge. The major drawback of this construction is its restriction 
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to spherically symmetric vector fields “attached” directly to the BH. We now show an 
alternative construction that allows for more generic field configurations. 


3.2. Condensates around black holes: Gaussian Initial Data solutions 

To access physically more interesting cases, such as non-spherically symmetric field 
profiles or “clouds” surrounding the BH in some distance, we construct initial data 
representing a (vector field) condensate around a BH, for Gaussian-like initial data 
(GID). After imposing (3.5) we still have the freedom to specify X^ and, as before, the 
momentum constraint (3.6) has the known Bowen-York solution |64j . Our strategy now 
is as follows: 


first we solve the “Gauss” constraint (3.8). We specify the “shape” of the condensate 
as 


X d 


= cTOe(^), 


(3.11) 


where 


R[r ) = exp 


r 0 y 


<j- 


Z(e,<p)=J2 C l ™Ylm(9,<p) 

Im 


composed of a Gaussian wavepacket R{r ) (of width a centered around r 0 ) with 
amplitude C. The angular distribution is given by a superposition of spherical 
harmonics Yi m (9, p) such that X^ is real. With this ansatz we can solve the resulting 
Poisson equation analytically for V and, hence, for the conformal field E l . 

Now we insert the ansatz for X^ and the corresponding V into the Hamiltonian 
constraint. Since we want to consider BH solutions we use the puncture 
approach [52J and set 


ip = VYl + u , 


( , M 

^BL — 1 + X - i 

2 r 


(3.12) 


where ip bl denotes the Brill-Lindquist conformal factor for a BH of mass M 
located at the origin and the correction u is a smooth function. The Hamiltonian 


constraint (3.7) becomes an elliptic equation for u which we solve numerically 
by modifying the TwoPUNCTURES code [67], which is part of the EINSTEIN 
Toolkit [68) [69]. Having been originally developed to calculate four-dimensional 


vacuum puncture data for binary BH configuration, TwoPUNCTURES has proven 
to work equally well when modified to handle different scenarios P31EQHZ2]. As we 
will see, this is also the case in our current configuration. 

We have to ensure that the source terms entering the Hamiltonian remain regular at the 
puncture and, furthermore, that the resulting elliptic PDE is well-posed. To convince 
ourselves of the latter property we linearize the equation around a known solution i/j 0 , 
ip — Vb + G with |e| <C ipo- Then, the linearized Hamiltonian is 

7 ,,, * 3 


77,i,, = Ae 


8ipo 


+ 


= Ae-he = 0. (3.13) 
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The existence of a unique solution requires h > omm which is clearly satisfied by 
our construction. 

To further simplify this procedure we start by considering time-symmetric 
configurations, i.e. A^j = 0. Note, that we are in no way limited to such configurations; 
we merely choose to focus first on these simpler cases before advancing to more involved 
examples. We will now explicitly write specific constructions that we have evolved. 


Gaussian initial data I — spherical symmetry: We first focus on spherically 
symmetric initial data and set Y(9,(p) = 2i/frloo = 1 and p — 1 in (3.11). With 
this ansatz we obtain 


** = ^f-R{r) , 

P r CqovW 

2 r 2 


a 


exp ( y ) — R(r ) + 


V^o 


cr 


(?) 


erf ( — ) + erf 


r -r 0 
a 


E e = 0 , E* = 0 , 


(3.14) 

,(3.15) 

(3.16) 


where R(r ) is given in (3.11). We obtain the “electric” field in Cartesian coordinates 
by a standard coordinate transformation. If we insert this solution into the source 
terms in (3.7) and recall that ^ ~ 1/r we can convince ourselves that they remain 
regular at the puncture. 

Gaussian initial data II — the dipole: We next construct initial data sourced by 
an axisymmetric dipole field, i.e., we set the angular function Y(9,p) = Y w (9,p) 
in (3.11) and choose the exponent p — 0. Specifically, we take the ansatz 

X* = C w R(r)Y w (6 , p) = C w J R(r) cos 9 , 

V 47t 

where R(r) has been defined in (3.11) and C w denotes the amplitude. We solve (3.8) 

a' . 1^-^—i 

with boundary conditions lirn^oo E l = 0 to find 


(3.17) 


E r = 


CloPyO- COsd | ^ ^ + r 2 + ^ _ 2(J ^2 + ^ exp 


a/487t r 

— y/n (r 3 + 2rg + 3r 0 a 2 ) erf 


a- 


r — r 0 
a 


a/vt r 0 (2 rl + 3cr 2 ) erf — + y/nr c 

2 a (r 2 + rr 0 + + a 2 ) R(r) — 2a (r^ + a 2 ) exp 


(3.18) 


p _ ^lo/iyo- sin 9 
Vl92v t r 4 

+y/n (2r 3 — 2rg — 3r 0 a 2 ) erf 
— y/nr 0 (2uq + 3cr 2 ) erf 

E v = 0 . 


a- 


r -r 0 
a 


ro 

l a 


2\/7rr 3 |, 


(3.19) 

(3.20) 


The source terms in the Hamiltonian constraint (3.7), given by our solution (3.17) 
and ( |3.18[ ), remain regular (and indeed vanish) at the puncture. We show in 
figure [I] an example of the correction term u (see equation (3.12)) obtained with 
this construction. 
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Figure 1. Typical example of the correction term u for a dipole configuration, 
equation (3.17), plotted along the x (left panel) and z (right panel) axis. It refers 
to a cloud with mass coupling M/x v = 0.4 centered around r$/M = 6 with width 
cr/M = 1 and amplitude C\q = 0.1. All dimensionful quantities are in units of the 
initial BH mass parameter M = 1. 


4. Wave extraction 

To obtain information about waves emitted in form of gravitational and electromagnetic 
radiation we employ the Newman-Penrose (NP) formalism [75]. We will use the notation 
conventions of [55]. In this spinor-inspired formalism the radiative degrees of freedom 
are given by a set of (complex) scalars. They are defined as contractions of the Weyl 
and Maxwell tensor, respectively, with the vectors of a null tetrad , #) with 

k^l^ = — 1 = We construct the null-tetrad vectors from 

e = 4= K + u") , 

y/2 

= —= ( v M + iw M ) , wE = -4= (u M — iwE , 
a/2 a/2 

where is the timelike unit normal vector and the spatial vectors (u l ,v\w l ) form a 
Cartesian orthonormal basis (see, e.g., (A2)-(A5) in [T6]). Asymptotically the basis 
vectors (u l ,v l ,w l ) behave as the unit radial, polar and azimuthal vectors. 

We extract the radiation carried in the vector fields by computing the NP scalars 

= W^rrf , = Uv iw (Pr + m^m") , <h 2 = , (4.1) 

which are reconstructed from the dynamical variables (A), E l ) using 

$o = - \ (E t v l + uV {DjXi - DiXjj) 

- l - (. EW + uV {D^ - DiXj)) , 

- ^v l w J (DjXi - DiXj) , 


k» = — OrE - uE , 

y/2 


(4.2) 

(4.3) 
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$2 = \ (Etf - uV {DjXi - DiXj)) 

- ( E % w i - Ew j (DjXi - DiXj)) . (4.4) 

The gravitational radiation is encoded in the in- and outgoing Weyl scalars (in an 
appropriate null tetrad) 

4'o = C pLVpa k p fh u k p rff , T 4 = C pvp J p m u l p rrf . (4.5) 

In analogy with electromagnetism it is convenient to decompose the Weyl tensor into 
its electric and magnetic components 


E- ■ = ry^-ryp C n v n° 

ij 1 i 1 j w per' 1 ' 1 


B = 'WV- 

til] 


• /"~Y V (7 


(4.6) 


where *C pupa = \e a ^ pa C pva p denotes the dual Weyl tensor. In terms of the evolution 
variables they are given by 

Bij = e(i\ kl D k A\j)i , (4.7) 

E^Rf-lA^.r+lKA, 

+ [EiEj - HyXiXj - 7 fcZ (D k Xi - DiX k ) (D^ - DjXi )] TF , (4.8) 


where is the trace-free part of the extrinsic curvature, [.] TF denotes the tracefree 
part (computed with the physical metric %j) and in the last relation we used (2.14). 
Finally, the Weyl scalars are computed from 


T 0 = 


T 4 = 


-ifjj (y l vi — + BijV l wi 

+ i ^ EijV l wi — - ( v l vi 

^Eij (v l iA — w l w — B, :J v 1 vj ] 

— i f EijV % vA + - ( v l vi — w % nP 


) 

) 



(4.9) 

(4.10) 


in terms of the electric and magnetic parts of the Weyl tensor in the Cartesian 
orthonormal basis (u l , w % ). 

At a given extraction radius R ex , we perform a multipolar decomposition by 
projecting T 4 , $ 4 and $2 onto spin-weighted spherical harmonics with s = —2, 0, 
and —1, respectively, 


*4(i, 8A) = ]T y m (i) - 2 VW(M) , (4.11) 

l,m 

*l(t, 0A)=Y, </>) , ( 4 - 12 ) 

l,m 


(4.13) 
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5. Numerical results 


5.1. Numerical setup and code description 

To simulate the interplay between BHs and massive vector fields we have implemented 
new thorns as part of the Lean code, originally presented in m for vacuum spacetimes. 
Lean uses the BSSN formulation of the Einstein equations [50j [5TJ with the moving 
puncture method J53] 54], is based on the Cactus Computational toolkit [75] . the 
Carpet mesh refinement package ESI (80] and uses AHFinderDirect for tracking 
apparent horizons [ST, 52]. We employ the method-of-lincs, where spatial derivatives 
are approximated by fourth-order finite difference stencils, and we use the 4 th -order 
Runge-Kutta scheme for the time integration. For further details on the numerical 
methods we refer the reader to mm 

The required extensions to Lean to accommodate for the Einstein-Proca system 
were coded independently by two of the authors (H.W. and M.Z.), thus allowing for 
independent cross-checks. Furthermore, since one recovers the standard Einstein- 


Maxwell case by setting ^y = 0 in the initial data construction outlined in section 3.1 


we 


have checked /xy = 0 evolutions against the corresponding (and tested) results obtained 
in [HI 183]. In all cases we have observed excellent agreement. 

We have evolved different initial profiles, as outlined in section |3] and monitored 
gauge invariant quantities such as the area of the BH’s apparent horizon (AH) and 
the NP scalars $ 0 , 1 , 2 , cf. (4.1). In table [l] we list the simulations done with the initial 


configuration outlined in section |3.1| where for fixed BH mass parameter M = 1 we can 
vary two parameters, C and /iy. Analogously, we list in table [2] simulations done with 
the initial configurations of section 332, where we have the parameters C 0 o, CT 0 , /xy, r 0 
and cr. 

To conclude the description of the implementations, we summarize their accuracy. 
We have simulated SID with a mass coupling M/xy = 0.4 and amplitude C/Ad = 0.4 
at three different resolutions h c /M = 1/84, h m /M = 1/92 and h^/M = 1/100. A 
convergence analysis reveals a numerical error in the waveforms 0?° of A(j)f° / (fi° h < 3% 
at early times (t/M = 200) which increases to A</>?°/0i°/i 1$ 9% after an evolution of 
t/M = 2000. 


5.2. Evolutions of Proca fields coupled to nonrotating BHs 

Our numerical results are summarized in figures [2]-[5] and the generic, qualitative features 
agree with previous studies on linearized Proca fields around BHs [43] and with nonlinear 
evolutions of massive scalars [34] , 

For example, figure [2] shows the evolution of SID data, with two different—and 
small—values of the amplitude C, for M /xy = 0.4. The two curves overlap after a linear 
re-scaling, showing that nonlinear effects are negligible for this data. More importantly, 
the time evolution generically shows the presence of a beating pattern. As explained 
previously [3U(43j, these are the result of the presence of multiple overtones with similar 
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Table 1. List of simulations performed with the SID of section |3.1[ along 
with parameters used. M always denotes the mass parameter in the initial data 
construction. The numerical grid structure used (in the notation of section II E of m 
was the following {(256,176,64,32,16,8,2,1,0.5), M/80}. Note, that we have actually 
run a larger set of simulations and present here only those that are shown later on. 


Run 

C/M 

M/iy 

Ei_c=0.4_mu=0.2 

0.4 

0.2 

Ei_c=0.05_mu=0.4 

0.05 

0.4 

Ei_c=0.4_mu=0.4 

0.4 

0.4 

Ei_c=0.4_rau=0.5 

0.4 

0.5 

Ei_c=0.4_mu=0.9 

0.4 

0.9 

Ei_c=0.4_rau=l.5 

0.4 

1.5 


Table 2. List of simulations performed with the GID of section |3.2[ along with 
parameters used. As before M = 1 is the initial mass parameter of the BH. The 
numerical grid structure used (in the notation of section II E of (77j) is given by 
{(256,176,64,32,16,8,2,1,0.5), M/80}. 

Run 


Coo/M 

Cio/M 

To/M 

a/M 

H\jM 

Madm/M 

phi_c00=0.1. 

_mu=0.2 

0.1 

0 

6.0 

1.0 

0.2 

1.00027 

phi_c00=l. 0. 

_mu=0.2 

1.0 

0 

6.0 

1.0 

0.2 

1.02668 

phi_c00=0.1. 

_mu=0.4_r0=5 

0.1 

0 

5.0 

1.0 

0.4 

1.00213 

phi c00=0.1. 

_mu=0.4_r0=12 

0.1 

0 

12.0 

1.0 

0.4 

1.00529 

phi_cl0=0.1. 

_mu=0.1 

0 

0.1 

6.0 

1.0 

0.1 

1.00011 

phi_clO=l. 0. 

_mu=0.1 

0 

1.0 

6.0 

1.0 

0.1 

1.01091 

phi_cl0=0.1. 

_mu=0.2 

0 

0.1 

6.0 

1.0 

0.2 

1.00051 

phi_cl0=0.1. 

_mu=0.4 

0 

0.1 

6.0 

1.0 

0.4 

1.00320 


frequencies (i.e., modes with different n in equation (2.26)). 

Similar results are obtained for other spherically symmetric configurations using 
SID data. In figure [3] we show the scalar profile </>5° for various mass couplings 
and keeping the amplitude C of the initial data fixed. In all cases we observe a 
(slowly) damped oscillatory pattern, corresponding to QBS-excitations (see section 2.4). 
Oscillation frequencies and damping times are here uniquely determined by the mass 
parameter /iy. For example, for M/iy = 1.5 the waveform is best fitted by a damped 
sinusoid of the form e'°' 021t sin 1.48t. This agrees very well with the fourth-overt one 
of Proca fields around a nonrotating BH, which we computed using frequency-domain 
calculations [HI HE SI] to be lu = 1.476 — A).02158. For smaller mass couplings, the 
signal is extremely long-lived, as can be seen from the top panel in figure [3] 

In fact, very long-lived, spherically symmetric configurations are one of the features 
borne out of our work. They arise for generic initial data. For example, GID data 
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M/x v =0.4 



0 200 400 600 800 1000 1200 

t/M 


Figure 2. cjy^ for the SID configurations of section 3.1 with M/zv = 0.4, extracted 
at i?ex = 140M. The (blue) solid line shows the case with amplitude C = 0.05, while 
the (green) dashed-dotted line refers to an amplitude C = 0.4. We have scaled the 
C = 0.05 results by 8 = 0.4/0.05; the two curves then overlap, confirming that our 
runs still he in the linear regime. The beating patterns in these curves are due to the 
presence of multiple overtones excited by the initial data. 


of section |3.2[ with more free parameters gives us a richer structure in the oscillation 
pattern, but also long-lived signals. This is more clearly seen in figure [4j where we 
plot (/>5° generated from SID and spherically symmetric GID initial configurations. 
The features worth mentioning are the beating pattern due to excitation of different 
overtones, but most specially the late time behaviour, which is independent of the 
initial data and is, basically, a very long-lived damped sinusoid. The signal can last for 
thousands of dynamical timescales, and the perturbative results indicate that long-lived 
modes might be present. 

In other words, QBS states are apparent in the evolution of the initial data. In 
some cases, it is possible to also observe late-time, oscillating power-law tails. While 
they have been buried under the long-lived QBS modes in the spherically symmetric 
cases, we do find the massive tails for dipolar initial configurations. We therefore focus 
now on simulations with type II GID data of section T2 In figure [5] we present a 


set of examples referring to the runs phi_cl0=0. l_mu=0.1, phi_cl0=0. l_mu=0.2 and 
phi_cl0=0. l_mu=0.4, where the dipole mode has been excited. Specifically, we display 
the l — l,m — 0 mode of the NP scalars </>]° and, for the smallest mass coupling, 
02° together with the envelopes for the power-law tails given in (2.28). The former 
set of waveforms is consistent with an intermediate-time, S = — 1,1 = 1 massive tail 
~ sin(/iyt). Furthermore, the waveform for the Mfi y = 0.4 case nicely exhibits the 
transition from the mode-dependent intermediate-time tail to the universal late-time 
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C=OAM, Mfi v =0.2 



0 500 1000 1500 2000 2500 3000 3500 4000 

(: t-RJ/M 


C=0AM, M^ v =0.4 





Figure 3. Results for SID configurations as described in section [3~T| with amplitude 
C = 0.4M and different mass-couplings, varying between M/j,y = 0.2 (top left) to 
Mfiy = 1.5 (bottom right). We plot the monopole waveform d? 0 , rescaled by the 
respective extraction radius. Due to spherical symmetry this is the only nontrivial 
mode excited throughout the evolution. In a massless /zy = 0 configuration, </>i° would 
measure the monopole component of the electric field, i.e., the electric charge. The 
inset in the last panel shows the data in a logarithmic scale with a fit of the form 
e -o. 02 i^ w j 1 j c ] rl matches very well with frequency-domain calculations. 


tail ~ £ -5 / 6 sin(/z v £), as can be seen in the bottom left panel of figure [HJ We expect a 
similar behaviour for the smaller mass-couplings for which this transition occurs after 
t 1000M, longer than our evolution time. Instead, the second NP scalar (pl° that we 
show exemplary for the small coupling case (see top left panel of figure [5]) corresponds 
to a different polarization for which the intermediate-time tail decays faster and, indeed, 
we observe an early onset of the universal late-time tail ~ t~ 5//6 sin(/iyt). 

To complete the picture illustrating the nonlinear interaction between massive 
vector fields and (nonrotating) BHs we have analysed the properties of the apparent 
horizon. We show the relative change in the AH area |Aah/Aah(£ = 0) — 11 in the 
bottom right panel of figure [5j This encodes the direct (nonlinear) response of the BH 
to the impinging held: In all models a fraction of the massive vector held is absorbed 
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Mfiy = 0.4 

- phi_c00=0.1_mu=0.4_r0=5 

- - phi_c00=0.1_mu=0.4_r0=12 

- - Ei c=0.4 mu=0.4 



1000 1500 2000 

t/M 


Figure 4. We present the monopole mode </q, rescaled by the extraction radius 
l?ex = 100M, for a Proca field with M/i y = 0.4. The different curves correspond to 
different initial configurations; specifically the (blue) solid and (green) short-dashed 


lines are type I GID (see section 3.2) where the spherically symmetric Gaussian was 
centered around r 0 = 5 M and ro = 12M, respectively, and the (red) long-dashed 
line refers to SID (see section 3.1). Although the BH’s immediate response varies 


substantially for different initial setups yielding more or less pronounced beating 
patterns, the overall oscillation frequency appears to be the same as can be seen in the 
inset. 


by the BH, the exact amount of which depends both on the initial energy content in 
the perturbing held and on its life-time. This can be seen from figure [5] by comparing 
the development of the models with small mass couplings (top panel) to the one with 
M/j ,v = 0.4 (bottom left panel). In the former cases the AH area increases upon the 
infall of the held after which it stays approximately constant consistent with the early 
onset of the power-law decay in the waveform. In the latter case the Proca held lingers 
around yielding a longer absorption phase and only when it reaches the tail-stage the 
AH area saturates to an approximately constant value. 

6. Final remarks 

In this work we have introduced a formalism to numerically integrate the Einstein-Proca 
system, describing General Relativity minimally coupled to a massive vector held, and 
we have used this formalism to begin to study how BHs evolve dynamically in this 
theory. Our results agree with perturbation theory expectations and pave the way for 
future nonlinear evolutions in more complex scenarios. Of particular interest, currently 
under development, is the study of configurations with rotating BHs. We find that 
the time-development of generic classes of initial data lead to extremely long-lived, 
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Figure 5. We present <j>\ 0 (solid lines) and $1° (dotted line in the top left 
panel) extracted at R ex = 100 M, obtained from an initial GID dipole configuration 
with M/i v = 0.1 (top left), M/xy = 0.2 (top right) and M/xy = 0.4 (bottom 
left). The (magenta) short-dashed and (red) long-dashed lines are envelopes of 
the (massive) power-law tails corresponding, respectively, to intermediate-time tails 
~ < -3 / 2 compatible with an S = —1,1 = 1 excitation, and to very late-time tails 
~ f -5 / 6 . In the bottom right panel we plot the relative change in the area of the 
corresponding AHs. The AH area increases during the first interaction between field 
and BH, and saturates to an approximately constant value when the power-law decay 
stage is reached. The downward trend observed at late times in the AH area is likely 
due to constraint violating modes, possibly due to reflections at the outer boundary of 
the computational domain (sitting at ~ 256 M). 

nearly periodic solutions, generalizing the findings of Ref. [S3, [S5| to the vector sector. 
Truly periodic, asymptotically flat hairy BHs have recently been ruled out |86j, but our 
results show that BHs in Proca theories might develop what for all purposes is a hair 
on timescales extremely large compared to any meaningful dynamical timescale. 
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Appendix A. BSSN equations 


For completeness, we display here the complete BSSN equations that we numerically 
integrate. 


{d t -Cp)X i 
(dt - Cp) Z 
(d t -Cp)X 4> 
(d t — Cp) E i 

(d t - Cp) la 
(d t - Cp) x 
(< d t -Cp)K 

(dt — Cp) Aij 
(d t - Cp) P 


= -ax - adiX^ - X^a 
= adiE 1 - ^ ax~ 1 E l d i x + - anZ 

= aKX, p - af j xdjXi + oyA,P + j d jX - X^Xdja - aZ 

= aKE i + afix^Xj + a X f j djZ + (d 3 X k - d k X 3 ) 

+ ax 2 f J X kl ( D k djXi - D k diX^j + (djXid k x - d k X 3 dix) 

2 aAtj , 

2 

= g a X K , 

= [...]+ 47 va(E + S ), 


= ...] — 877a 




= 167ray \f , 


(A.l) 

(A.2) 

(A.3) 

(A.4) 

(A.5) 

(A.6) 

(A.7) 

(A.8) 

(A.9) 



REFERENCES 


21 


where [...] denotes the standard right-hand side of the BSSN equations in the absence 
of source terms; the source terms are determined by 


E =n a n^T a p, (A. 10) 

3i =-li a ^T af} , (A. 11) 

(A. 12) 

5 = 7 VSij. (A. 13) 


It is convenient to introduce the “magnetic” field B l 
B i = e^ k DjX k = X 3/2 ef (djX, + \x~\XAx + X k d jX ~ 7 ron 7ifc*nAx)) , (A. 14) 

as short-hand, where ef 1 ' is a tensor density taking the values of (0, ±1) as in flat space. 


We then have 

4t t{E + S) = x~ X la {B i E j + B'Bi) + 2&X} , (A.15) 

4tt (xSa ~ f 7y) = - lik%iX~ X {E k E l + B k B l ) + 1 -f{E 2 + B 2 ) + /4xAX, 

-$r^xfX k X u (A.16) 

4 v X ~ X f = 7 li X~ 3/2 ef jk EiB k + ^X^fX k , (A.17) 

87 rE = x~%3 {E i E j + B l B>) + i4 (A| + X l kl X k X x ) . (A. 18) 
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